Chapter 2.12 - Circular Restricted Three-Body Problem¶

In this section, we solve the three-body problem, subject to some restrictions. In particular:

There are two primary masses, and the mass of the tertiary object is extremely small in comparison to \(m_1\) and \(m_2\)
The mass of \(m_1\) is larger than \(m_2\)
The two primary objects orbit in a circle around their center of mass

Although these assumptions seem fairly restrictive, they actually represent several very important physical situations: the Earth-Moon system, as well as the orbits of many of the planets around the sun, with a man-made object as the third mass! This is called the Circular Restricted Three-Body Problem (CRTBP or CR3BP), because the orbits are restricted to circles and the mass of the third body is restricted to be much smaller than the other two.

The orbit of the moon around the earth is approximately circular, with a mean eccentricity of 0.054, and semimajor and semiminor axes of 384,400 km and 383,800 km, respectively. The center of mass of the system occurs at a distance of 4,600 km from the Earth’s center, about 72% of the radius of the Earth. This data comes from Wikipedia.

Similarly, the orbits of Venus, Earth, Jupiter, Saturn, Neptune, and Uranus around the sun all have eccentricities less than 0.06, according to the NASA Planetary Fact Sheet.

Unlike the two-body problem, there is no general closed-form solution to this problem. Closed-form means an analytical equation we can solve. But we can construct the equations of motion to find some interesting parameters of the orbits.

Orbit of Primary Masses¶

We first attach a non-inertial coordinate system to the barycenter of the system of \(m_1\) and \(m_2\), such that the \(x\)-axis of this coordinate system points towards \(m_2\). The distance from \(m_1\) to \(m_2\) is \(r_{12}\), which is also the radius of the circular orbit.

The \(y\)-axis of this coordinate system is in the orbital plane, and the \(z\)-axis is perpendicular to the orbital plane, in the same direction as the angular momentum vector. In the rotating reference frame, \(m_1\) and \(m_2\) appear to be stationary.

The inertial angular velocity is given by

\[\vector{\Omega} = \Omega\uvec{k}\]

where

\[\Omega = \frac{2\pi}{T}\]

and the orbital period is for a circular orbit:

\[T = \frac{2\pi}{\sqrt{\mu}}r_{12}^{3/2}\]

Plugging this in for \(\Omega\), we find:

\[\Omega = \sqrt{\frac{\mu}{r_{12}^3}}\]

where the gravitational parameter is given by:

\[\mu = GM = G\left(m_1 + m_2\right)\]

Next, we want to find the positions of the two masses relative to the barycenter. By definition, the two masses lie in the orbital plane, so their \(z\)-coordinate is going to be zero. Since the line that connects \(m_1\) and \(m_2\) goes through the barycenter, then their \(y\)-coordinates must be zero as well.

In that case, we only need to find the \(x\)-coordinates, which we can do from the equation for the center of mass:

\[m_1 x_1 + m_2 x_2 = 0\]

We need a second equation to solve for \(x_1\) and \(x_2\). Since \(m_1\) is to the left of the center of mass, \(x_1\) is negative. Then \(x_2\) is the distance from \(m_1\) to \(m_2\) plus the distance from \(m_1\) to the barycenter:

\[x_2 = x_1 + r_{12}\]

Then we define two dimensionless ratios:

\[\begin{aligned}\pi_1 &= \frac{m_1}{m_1 + m_2} & \pi_2 &= \frac{m_2}{m_1 + m_2}\end{aligned}\]

and solve for \(x_1\) and \(x_2\) in terms of these:

\[\begin{aligned}x_1 &= -\pi_2 r_{12} & x_2 &= \pi_1 r_{12}\end{aligned}\]

Orbit of the Tertiary Mass¶

Now let’s add the much smaller, tertiary mass into the system. We’ll use the symbol \(m\) for this mass, without a subscript. We want the equation of motion, that is, Newton’s second law. For that we need the acceleration, which we will derive from the position.

Position, Velocity, and Acceleration¶

The position of the tertiary mass relative to the barycenter is:

\[\vector{r} = x\uvec{\imath} + y\uvec{\jmath} + z\uvec{k}\]

The position of the tertiary mass relative to \(m_1\) is:

\[\vector{r}_1 = \left(x - x_1\right)\uvec{\imath} + y\uvec{\jmath} + z\uvec{k} = \left(x + \pi_2 r_{12}\right)\uvec{\imath} + y\uvec{\jmath} + z\uvec{k}\]

and finally, the position of \(m\) relative to \(m_2\) is:

\[\vector{r}_2 = \left(x - \pi_1 r_{12}\right)\uvec{\imath} + y\uvec{\jmath} + z\uvec{k}\]

Then we want the inertial velocity of \(m\), so we need to account for the rotating frame of reference, so:

\[\dot{\vector{r}} = \vector{v}_G + \vector{\Omega}\cross\vector{r} + \vector{v}_{\text{rel}}\]

where \(\vector{v}_G\) is the absolute velocity of the barycenter and \(\vector{v}_{\text{rel}}\) is the velocity calculated in the moving coordinate system:

\[\vector{v}_{\text{rel}} = \dot{x}\uvec{\imath} + \dot{y}\uvec{\jmath} + \dot{z}\uvec{k}\]

Then we can find the absolute acceleration of \(m\):

(4)¶\[\ddot{\vector{r}} = \vector{a}_G + \dot{\vector{\Omega}}\cross\vector{r} + \vector{\Omega}\cross\left(\vector{\Omega}\cross\vector{r}\right) + 2\vector{\Omega}\cross\vector{v}_{\text{rel}} + \vector{a}_{\text{rel}}\]

This equation can be simplified because we showed that the acceleration of the barycenter is zero for the two-body problem, \(\vector{a}_G = 0\). In addition, the angular velocity is constant since the orbit is circular, so \(\dot{\vector{\Omega}} = 0\). Then, Eq. (4) can be simplified to:

\[\ddot{\vector{r}} = \vector{\Omega}\cross\left(\vector{\Omega}\cross\vector{r}\right) + 2\vector{\Omega}\cross\vector{v}_{\text{rel}} + \vector{a}_{\text{rel}}\]

where

\[\vector{a}_{\text{rel}} = \ddot{x}\uvec{\imath} + \ddot{y}\uvec{\jmath} + \ddot{z}\uvec{k}\]

Plugging everything in and simplifying:

\[\ddot{\vector{r}} = \left(\ddot{x} - 2\Omega\dot{y} - \Omega^2 x\right)\uvec{\imath} + \left(\ddot{y} + 2\Omega\dot{x} - \Omega^2 y\right)\uvec{\jmath} + \ddot{z}\uvec{k}\]

Newton’s Second Law¶

For the tertiary body, the forces are due to both of the other masses:

\[m \ddot{\vector{r}} = \vector{F}_1 + \vector{F}_2\]

where \(\vector{F}_1\) is the force from \(m_1\) on \(m\) and \(\vector{F}_2\) is the force from \(m_2\) on \(m\).

The two forces are found by Newton’s law of gravitation:

\[\begin{split}\begin{aligned}\vector{F}_1 &= -G\frac{m_1 m}{r_1^2} \left.\uvec{u}_r\right)_1 = -\frac{\mu_1 m}{r_1^3}\vector{r}_1 \\ \vector{F}_2 &= -G\frac{m_2 m}{r_2^2}\left.\uvec{u}_r\right)_2 = -\frac{\mu_2 m}{r_2^3}\vector{r}_2\end{aligned}\end{split}\]

where

\[\begin{aligned}\mu_1 &= G m_1 & \mu_2 &= G m_2\end{aligned}\]

and

\[\begin{aligned}\left.\uvec{u}_r\right)_1 &= \frac{\vector{r}_1}{r_1} & \left.\uvec{u}_r\right)_2 &= \frac{\vector{r}_2}{r_2}\end{aligned}\]

Dividing through by \(m\), we find:

\[\ddot{\vector{r}} = -\frac{\mu_1}{r_1^3}\vector{r}_1 - \frac{\mu_2}{r_2^3}\vector{r}_2\]

Now we substitute for \(\ddot{\vector{r}}\) and split out by components to have three scalar equations of motion for the CRTBP:

\[\begin{split}\begin{aligned}\ddot{x} - 2\Omega\dot{y} - \Omega^2 x &= -\frac{\mu_1}{r_1^3}\left(x + \pi_2 r_{12}\right) - \frac{\mu_2}{r_2^3}\left(x - \pi_1 r_{12}\right) \\ \ddot{y} + 2\Omega\dot{x} - \Omega^2 y &= -\frac{\mu_1}{r_1^3}y - \frac{\mu_2}{r_2^3}y \\ \ddot{z} &= -\frac{\mu_1}{r_1^3}z - \frac{\mu_2}{r_2^3}z\end{aligned}\end{split}\]

There are a few things we can note from these equations. First, the \(x\) and \(y\) equations are coupled; the \(x\) equation depends on \(y\), and the \(y\) equation depends on \(x\). However, the \(z\) equation is decoupled from the other two, so if \(m\) starts in planar motion, it will remain there.

Nondimensional Equations of Motion¶

Next, let’s make these equations nondimensional. This offers the advantage of being general for any system we want to study, and removing the dependence on the rate of rotation of the coordinate system. We start by defining a few characteristic parameters.

First, the characteristic distance will be \(r_{12}\) and the characteristic mass will be \(\pi_2\), for convenience. We could choose \(\pi_1\) without loss of generality. We will also define a characteristic time:

\[t^* = \sqrt{\frac{r_{12}^3}{G\left(m_1 + m_2\right)}}\]

This characteristic time is used to generate the nondimensional time:

\[\tau = \frac{t}{t^*}\]

Then we need to define the dimensionless position vectors. The vector from \(G\) becomes:

\[\vector{\rho} = \frac{\vector{r}}{r_{12}} = x^*\uvec{\imath} + y^*\uvec{\jmath} + z^*\uvec{k}\]

where \(x^* = x/r_{12}\), and similar for \(y^*\) and \(z^*\).

The other two position vectors, from \(m_1\) and \(m_2\), respectively, become:

\[\begin{split}\begin{aligned}\vector{\sigma} &= \frac{\vector{r}_1}{r_{12}} = \left(x^* + \pi_2\right)\uvec{\imath} + y^*\uvec{\jmath} + z^*\uvec{k} \\\vector{\psi} &= \frac{\vector{r}_2}{r_{12}} = \left(x^* - 1 + \pi_2\right)\uvec{\imath} + y^*\uvec{\jmath} + z^*\uvec{k}\end{aligned}\end{split}\]

Note that \(\pi_1 = 1 - \pi_2\).

Making the acceleration equation nondimensional involves replacing both the \(\vector{r}\) and the \(t\) terms in the \(d^2\vector{r}/dt^2\):

\[\ddot{\vector{\rho}} = \frac{d^2\vector{r}}{dt^2}\frac{\left(t^{*}\right)^2}{r_{12}} = \frac{d^2\vector{\rho}}{d\tau^2}\]

Making the terms on the right hand side of Eq. (4) nondimensional is also the result of multiplying by \(\left(t^*\right)^2/r_{12}\). Note that the dimensions of \(\Omega\) are \(t^{-1}\):

\[\ddot{\vector{\rho}} = \left(\ddot{x}^* - 2\dot{y}^* - x^*\right)\uvec{\imath} + \left(\ddot{y}^* + 2\dot{x}^* - y^*\right)\uvec{\jmath} + \ddot{z}^*\uvec{k}\]

Now we have the nondimensional inertial acceleration, and we need to make Newton’s second law nondimensional:

\[\ddot{\vector{\rho}} = -\frac{1 - \pi_2}{\sigma^3}\vector{\sigma} - \frac{\pi_2}{\psi^3}\vector{\psi}\]

where \(\sigma = \mag{\vector{\sigma}}\) and \(\psi = \mag{\vector{\psi}}\). Now we can break this into the three scalar components as before:

\[\begin{split}\begin{aligned}\ddot{x}^* - 2\dot{y}^* - x^* &= -\frac{1 - \pi_2}{\sigma^3}\left(x^* + \pi_2\right) - \frac{\pi_2}{\psi^3}\left(x^* - 1 + \pi_2\right) \\ \ddot{y}^* + 2\dot{x}^* - y^* &= -\frac{1 - \pi_2}{\sigma^3} y^* - \frac{\pi_2}{\psi^3}y^* \\ \ddot{z}^* &= -\frac{1 - \pi_2}{\sigma^3}z^* - \frac{\pi_2}{\psi^3}z^*\end{aligned}\end{split}\]

There are three main advantages of this formulation of the equations of motion:

The equations are independent of the two masses \(m_1\) and \(m_2\), depending only on their relative sizes via \(\pi_2\)
The equations are independent of the rate of rotation of the reference frame
The equations are independent of the separation distance between \(m_1\) and \(m_2\), so can be used to represent any system

Lagrange Points¶

With the equations of motion, we can determine solutions by integrating them in time. As we will see, this can be a little bit tricky, depending on the situation. We also cannot solve the equations analytically, either in the dimensional or nondimensional case.

However, there are five equilibrium points in this system of equations that we can find. An object placed at one of the equilibrium points will, if it is not perturbed, remain at that location indefinitely. These locations are called Lagrange Points, although Lagrange himself called them libration points.

By definition, the acceleration and velocity components at the Lagrange points are zero:

\[\dot{x} = \dot{y} = \dot{z} = 0 \qquad \text{and} \qquad \ddot{x} = \ddot{y} = \ddot{z} = 0\]

and likewise for the nondimensional coordinates as well. Plugging these conditions in to the equations of motion, we find (taking the non-dimensional case):

\[\begin{split}\begin{aligned}-x^* &= -\frac{1 - \pi_2}{\sigma^3}\left(x^* + \pi_2\right) - \frac{\pi_2}{\psi^3}\left(x^* - 1 + \pi_2\right) \\-y^* &= -\frac{1 - \pi_2}{\sigma^3} y^* - \frac{\pi_2}{\psi^3}y^* \\ 0 &= \left(-\frac{1 - \pi_2}{\sigma^3} - \frac{\pi_2}{\psi^3}\right)z^*\end{aligned}\end{split}\]

In the last equation, for \(z^*\), the terms in the parentheses are both greater than zero. Therefore, there is no way for them to cancel each other out, and the only way for the equation to be satisfied is if \(z^* = 0\). Thus, the Lagrange points lie in the orbital plane.

Now, we need to solve for the \(x^*\) and \(y^*\) coordinates of the Lagrange points. There are two conditions that we need to consider, since they end up with different solutions:

\(y^* \neq 0\), giving the so-called equilateral Lagrange points
\(y^* = 0\), giving the so-called collinear Lagrange points

Equilateral Lagrange Points¶

To find the equilateral Lagrange points, we assume the \(y^* \neq 0\). Now we have 2 equations and 2 unknowns. The \(y^*\) in the second equation of motion will cancel from both sides, and we end up with:

\[1 = \frac{1 - \pi_2}{\sigma_3} + \frac{\pi_2}{\psi^3}\]

In the equation for \(x^*\) we see the \(\left(1 - \pi_2\right)/\sigma^3\) term, so let’s solve the \(y^*\) equation for that:

\[\frac{1 - \pi_2}{\sigma^3} = 1 - \frac{\pi_2}{\psi^3}\]

and plugging this into the \(x^*\) equation, we find:

\[x^* = \left(1 - \frac{\pi_2}{\psi^3}\right)\left(x^* + \pi_2\right) + \frac{\pi_2}{\psi^3}\left(x^* - 1 + \pi_2\right)\]

Simplifying, we find:

\[\psi^3 = 1\]

or, plugging in the definition of \(\psi\):

\[\psi = \frac{\mag{\vector{r}_2}}{r_{12}} \Rightarrow r_2 = r_{12}\]

Plugging the result of \(\psi^3 = 1\) back into the simplified \(y^*\) equation, we find:

\[\sigma^3 = 1\]

or, plugging in the definition of \(\sigma\):

\[\sigma = \frac{\mag{\vector{r}_1}}{r_{12}} \Rightarrow r_1 = r_{12}\]

Since if \(a = b\) and \(b = c\), then \(a = c\), we find:

\[r_1 = r_2 = r_{12}\qquad \text{and}\qquad \psi = \sigma\]

for the equilateral Lagrange points. Thus, the distance from \(m_1\) to \(m\) is the same as the distance from \(m_2\) to \(m\) and from \(m_1\) to \(m_2\). This defines an equilateral triangle, giving these Lagrange points their name.

From the definition of \(\vector{\sigma}\) and \(\vector{\psi}\), we can take their magnitudes and equate them to solve for the value of \(x^*\) at the equilibrium points:

\[\begin{split}\begin{aligned}\sigma^2 &= \mag{\vector{\sigma}}^2 = \left(x^* + \pi_2\right)^2 + \left(y^*\right)^2\\\psi^2 &= \mag{\vector{\psi}}^2 = \left(x^* - 1 + \pi_2\right)^2 + \left(y^*\right)^2\end{aligned}\end{split}\]

where we have also used the fact that \(z^* = 0\) for these points. We find that:

\[x^* = \frac{1}{2} - \pi_2\]

Plugging this back into the equation for \(\sigma^2\), we find:

\[y^* = \pm \frac{\sqrt{3}}{2}\]

The equilateral Lagrange points are given the symbols \(L_4\) and \(L_5\), for the positive and negative \(y^*\) values, respectively:

\[\begin{split}\begin{aligned}&L_4: & x^* &= \frac{1}{2} - \pi_2 & y^* &= \frac{\sqrt{3}}{2}\\&L_5: & x^* &= \frac{1}{2} - \pi_2 & y^* &= -\frac{\sqrt{3}}{2}\end{aligned}\end{split}\]

Collinear Lagrange Points¶

For the collinear Lagrange points, we set \(y = z = y^* = z^* = 0\) in our equations of motion. We make this choice by inspection, seeing that setting \(y = y^* = 0\) is one possible solution (the other case, \(y\neq 0\) we just handled).

The first equation of motion is then:

\[x^* = \frac{1 - \pi_2}{\sigma^3}\left(x^* + \pi_2\right) + \frac{\pi_2}{\psi^3}\left(x^* - 1 + \pi_2\right)\]

where, with \(y^* = 0\), \(\sigma^3\) and \(\psi^3\) are:

\[\begin{aligned}\sigma^3 &= \left\lvert x^* + \pi_2\right\rvert ^3 & \psi^3 &= \left\lvert x^* - 1 + \pi_2\right\rvert ^3\end{aligned}\]

Substituting these into the equation of motion:

\[0 = x^* - \frac{1 - \pi_2}{\left\lvert x^* + \pi_2\right\rvert ^3}\left(x^* + \pi_2\right) - \frac{\pi_2}{\left\lvert x^* - 1 + \pi_2\right\rvert ^3}\left(x^* - 1 + \pi_2\right)\]

We cannot cancel any of the terms because we don’t know the sign of the cubic term on the bottom, since we took the square root of the respective vectors to find the magnitude. In any case, this equation is cubic in \(x^*\), so it will have three roots for the nontrivial cases where \(0 < \pi_2 < 1\). In the cases where \(\pi_2 = 0\) or \(\pi_2 = 1\), either \(m_2\) or \(m_1\) must be zero, which aren’t very interesting.

There is no analytical solution to this equation, it must be solved numerically. There are many methods to solve the equation numerically, my suggestions are to use scipy.optimize.newton() in Python and fzero in Matlab. We will see how to use these functions in the next example.

In the meantime, it is useful to plot the solutions of this equation:

../_images/12-circular-restricted-three-body-problem_1_0.png

On this figure, the \(x\) axis is \(\pi_2\) and the \(y\) axis is \(x^*\). Remember that the nondimensional position of \(x^*_{m_1} = -\pi_2\) and of \(x^*_{m_2} = 1 + \pi_2\). For a given value of \(\pi_2\), the mass ratio of \(m_2\) to \(m_1 + m_2\), we can locate the positions of \(m_1\) and \(m_2\) from the dashed gray lines on the graph.

There is also an S-curve shape which represents the solution of the equation for the \(x^*\) position of the collinear Lagrange points. For a given value of \(\pi_2\), we can see there are 3 solutions of the function, corresponding to the three Lagrange points for that system.

By convention, the Lagrange points are numbered such that \(L_1\) lies between \(m_1\) and \(m_2\), \(L_2\) lies to the right of \(m_2\), and \(L_3\) lies to the left of \(m_1\). Thus, we can see on the figure that the upper part of the S-curve is the solution for \(L_2\). Below \(x^{*} = 1.0\), the solution is for \(L_1\), since that lies between \(m_1\) and \(m_2\). Finally, below \(x^{*} = -1.0\), the solution is for \(L_3\).

The figure below plots the five Lagrange points in nondimensional coordinates as a function of the mass ratio \(\pi_2\):

Lagrange Point Stability¶

As we discussed previously, gravity is a conservative force. As such, the force field can be defined in terms of a potential energy function. For the CR3BP, the pseudo-potential function in the rotating frame is given by [KLMR11]:

\[U(x^*, y^*) = -\frac{1 - \pi_2}{\sigma} - \frac{\pi_2}{\psi} - \frac{1}{2}\left[\left(1 - \pi_2\right)\sigma^2 + \pi_2 \psi^2\right]\]

A plot of this function is shown below, including the positions of the five Lagrange points, for \(\pi_2 = 0.3\):

../_images/12-circular-restricted-three-body-problem_6_0.png

Using this figure, we can get a qualitative sense of the stability of the Lagrange points. Imagine that we turn the potential function upside-down, and put a marble on each of the Lagrange points. We can see that \(L_4\) and \(L_5\) are at the bottom of a bowl. Slight displacements of the marble will cause it to return to the initial position, so these points are considered stable Lagrange points.

Since the equilateral Lagrange points are stable, objects placed in a small orbit, usually called a halo orbit, around those points will tend to remain there. The criterion for stability is:

\[\frac{m_1}{m_2} + \frac{m_2}{m_1} \geq 25\]

which will be satisfied if \(m_1/m_2>24.95994\) or \(\pi_2 < 0.0385209\). In the Earth-Moon system, that ratio is \(m_1/m_2 \approx 81.3\). However, in the Earth-Moon system, the \(L_4\) and \(L_5\) points are slightly destabilized by the influence of the sun. Nonetheless, there are clouds of dust which have collected at these points.

However, other pairs of \(m_1\) and \(m_2\) do have somewhat more stable \(L_4\) and \(L_5\) points, in particular, the orbit of Jupiter around the sun. There are groups of asteroids, called trojan asteroids that cluster around the stable Lagrange points in the orbit of Jupiter.

../_images/12-circular-restricted-three-body-problem_8_0.png

As you can see, there are clusters of asteroids at angles of \(\pm 60^{\circ}\) from Jupiter, corresponding to the \(L_4\) and \(L_5\) Lagrange points!

On the other hand, \(L_1\), \(L_2\), and \(L_3\) are all saddle points, meaning that the function increases when going in one axis, but decreases going in the other axis. This means that the three collinear Lagrange points are unstable and an object placed at one of those points, if perturbed, will diverge from the position.

Nonetheless, these are quite useful points for observation of the solar system. Several satellites have been placed at the \(L_1\) point of the Earth-Sun system for solar observation, and the James Webb Space Telescope is planned to launch to the \(L_2\) of the Earth-Sun system sometime this year.

\(L_1\) and \(L_2\) in the Earth-Sun system are about 1.5 million km towards the Sun and away from the Sun, starting at the Earth, respectively. \(L_3\) lies on the other side of the Sun, and has long been the predicted location of a hidden planet, since it could not be observed from Earth prior to the advent of satellite observation. Now, of course, we know there is no planet at that location.

The Jacobi Constant¶

Let’s return now to the graph of the potential function. The potential function represents the potential energy that the tertiary mass will have if it is located at a given \(x\)-\(y\) point in the orbital plane. The tertiary mass will also have some velocity, \(v\), with a corresponding kinetic energy.

From the law of conservation of energy, we know that the sum of the kinetic and potential energies of the tertiary mass will be constant. The potential energy (repeated here from above) is:

\[U(x^*, y^*) = -\frac{1 - \pi_2}{\sigma} - \frac{\pi_2}{\psi} - \frac{1}{2}\left[\left(1 - \pi_2\right)\sigma^2 + \pi_2 \psi^2\right]\]

The first two terms in this equation are the gravitational potential energy due to the position of the tertiary mass relative to \(m_1\) and \(m_2\). The third term is the potential energy of the centrifugal force induced by the rotation of the coordinate system.

In nondimensional coordinates, the mass-specific kinetic energy is:

\[\frac{1}{2} \left(v^*\right)^2 = \frac{1}{2}\left[\left(\dot{x}^*\right)^2 + \left(\dot{y}^*\right)^2\right]\]

Combining these equations, we find from conservation of energy:

(5)¶\[\frac{1}{2} \left(v^*\right)^2 - \frac{1 - \pi_2}{\sigma} - \frac{\pi_2}{\psi} - \frac{1}{2}\left[\left(1 - \pi_2\right)\sigma^2 + \pi_2 \psi^2\right] = C\]

The constant \(C\) is called the Jacobi Constant, and represents the total energy of the tertiary mass relative to the rotating reference frame.

Since \(C\) is a constant, the total energy of the tertiary mass is fixed. Consider the tertiary mass at some location \(\left(x^*_1, y^*_1)\) such that it has potential energy \(U_1\). Assume first that the velocity is zero. Then, the Jacobi constant is just equal to \(U_1\), and the mass cannot “climb” any higher out of the potential energy surface. Thus, there is a region of space where the mass cannot access because it doesn’t have enough energy!

Assume now that the velocity of the mass is \(v^*_1\). Then the Jacobi constant is equal to the sum of the kinetic and potential energies. If the mass wants to climb up the potential energy surface, it can trade kinetic energy for potential energy. Eventually, however, the kinetic energy and the velocity will go to zero, and the mass cannot climb any higher!

Now, let’s turn this problem around. We want to know, for a given value of \(C\), what regions of space will be inaccessible. Consider the tertiary mass with a certain value of \(C\) and at a particular location. As the mass moves, it exchanges energy between kinetic energy (velocity) and potential energy. At some position, the \(C\) will be equal to \(U\), and the velocity will be (by definition) zero. Thus, the mass cannot travel any further in that direction.

For a given value of the Jacobi constant, we can calculate the contours of zero velocity positions by setting \(v^* = 0\) in Eq. (5):

\[\frac{2 \left(1 - \pi_2\right)}{\sigma} + \frac{2 \pi_2}{\psi} + \left[\left(1 - \pi_2\right)\sigma^2 + \pi_2 \psi^2\right] + 2C = 0\]

Since the first three terms on the left are all positive, zero velocity curves correspond to negative values of the Jacobi constant. The figure below plots the forbidden regions, shown as shaded areas, for several values of \(C\):

../_images/12-circular-restricted-three-body-problem_10_0.png

Notes on Orbital Mechanics for Engineering Students

Chapter 2.12 - Circular Restricted Three-Body Problem¶

Orbit of Primary Masses¶

Orbit of the Tertiary Mass¶

Position, Velocity, and Acceleration¶

Newton’s Second Law¶

Nondimensional Equations of Motion¶

Lagrange Points¶

Equilateral Lagrange Points¶

Collinear Lagrange Points¶

Lagrange Point Stability¶

The Jacobi Constant¶